
package scalaSci

import scalaSci.math.LinearAlgebra.LinearAlgebra
import scalaSci.math.array.DoubleArray
import java.util._
import Jama._

object DoubleDoubleArr {


// construct an Array[Array[Double]]from a String 
// var dd = MDD("3.4 -6.7; -1.2 5.6")
def MDD(s: String): Array[Array[Double]] =  {
    var nRows = 1
    var nCols = 0
    for (i<-0 to s.length-1)   // count how many rows are specified
      if  (s(i)==';')
        nRows += 1

  // seperate rows to an ArrayBuffer
   var buf = new scala.collection.mutable.ArrayBuffer[String]()
   var strtok = new java.util.StringTokenizer(s, ";")  // rows are separated by ';'
   while (strtok.hasMoreTokens) {
         val tok = strtok.nextToken
         buf +=tok
      }    

// count how many numbers each row has. Assuming that each line has the same number of elements
 val firstLine = buf(0)
 strtok = new java.util.StringTokenizer(firstLine, ", ")  // elements are separated by ',' or ' '
 while (strtok.hasMoreTokens) {
   val tok = strtok.nextToken
   nCols += 1  
}

// allocate double array
var elems = new Array[Array[Double]](nRows, nCols)
  // pass elements from ArrayBuffer to double array
for (k <- 0 to buf.length-1)  {  // read array
   var currentLine = buf(k)
   strtok = new java.util.StringTokenizer(currentLine, ", ")  // elements are separated by ',' or ' '
var c=0 
while (strtok.hasMoreTokens) {  // read row
   val tok = strtok.nextToken
   var cNumber = tok.toDouble
   elems(k)(c) = cNumber
   c += 1
     }   // read row
   }  // read array

 // return the converted Array[Array[Double]]
 elems
 }  


 def reshape(darr: Array[Array[Double]], n: Int, m: Int): Array[Array[Double]] =  {
    var x_dy = darr.length
    var x_dx = darr(0).length
    if  ((x_dy*x_dx) != (n*m))  { return null }
    // create Matrix
    var values = new Array[Array[Double]](n, m)
    var yii = 0
    var xii = 0
    for (xi <- 0  to  m-1)
       for (yi <- 0  to  n-1)
          values(yi)(xi) = darr(yii)(xii)
      yii += 1
      if (yii >= x_dy)  {
         yii = 0
         xii += 1
      }
      return  values
     
 }


// columnwise min
def min(darr: Array[Array[Double]]): Array[Double] = {
    var N = darr.length;     var M = darr(0).length
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     var mn = darr(0)(ccol)  
     for (crow<- 1 to N-1)
        {
       var tmp = darr(crow)(ccol)
       if (tmp < mn)  mn = tmp
       }
     res(ccol) = mn
     }
    res
}

// columnwise max
def max(darr: Array[Array[Double]]): Array[Double] = {
    var N = darr.length;     var M = darr(0).length
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     var mx = darr(0)(ccol)  
     for (crow<- 1 to N-1)
        {
       var tmp = darr(crow)(ccol)
       if (tmp > mx)  mx = tmp
       }
     res(ccol) = mx
     }
    res
}

// columnwise sum
def sum(darr: Array[Array[Double]]): Array[Double] = {
    var N = darr.length;     var M = darr(0).length
    var sm = 0.0
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     sm=0.0
     for (crow<- 0 to N-1)
       sm += darr(crow)(ccol)
     res(ccol) = sm 
     }
    res
}
 
// columnwise mean
def mean(darr: Array[Array[Double]]): Array[Double] = {
    var N = darr.length;     var M = darr(0).length
    var sm = 0.0
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     sm=0.0
     for (crow<- 0 to N-1)
       sm += darr(crow)(ccol)
     res(ccol) = sm/N
     }
    res
}

// columnwise product
def prod(darr: Array[Array[Double]]): Array[Double] = {
    var N = darr.length;     var M = darr(0).length
    var pd = 1.0
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     pd=1.0
     for (crow<- 0 to N-1)
       pd *= darr(crow)(ccol)
     res(ccol) = pd 
     }
    res
}

 def Fill(n: Int, m:Int, vl: Double): Array[Array[Double]] = {
       var  v = new Array[Array[Double]](n, m)
       var i=0; var j=0;
       while  (i< n) {
            j=0
            while (j < m) {
              v(i)(j) = vl
              j += 1
              }
            i += 1
         }
          v
    }



def  ones(n:Int, m:Int): Array[Array[Double]] = {
     var  v = new Array[Array[Double]](n, m)
       var i=0; var j=0;
       while  (i< n) {
            j=0
            while (j < m) {
              v(i)(j) = 1.0
              j += 1
              }
            i += 1
         }
          v
    }

def  zeros(n:Int, m:Int): Array[Array[Double]] = {
     var  v = new Array[Array[Double]](n)
       var c=0
       while (c<n) {
          v(c) = new Array[Double](m)
          c += 1
       }
       var i=0; var j=0;
       while  (i< n) {
            j=0
            while (j < m) {
              v(i)(j) = .0
              j += 1
              }
            i += 1
         }
          v
    }


def rand(n: Int, m:Int): Array[Array[Double]] = {
      var  v =new Array[Array[Double]](n, m)
       var i=0; var j=0;
       while  (i< n) {
           j=0
           while (j < m) {
              v(i)(j) = java.lang.Math.random
              j += 1
              }
            i += 1
         }
        v
    }

def  eye(n:Int, m:Int): Array[Array[Double]] = {
     var  v = new Array[Array[Double]](n, m)
     var minCoord = n
     if ( m<minCoord)
        minCoord = m
     for (i<- 0 to minCoord-1)
         v(i)(i) = 1.0
     v
    }

def  eye(n:Int): Array[Array[Double]] =  {
    eye(n, n)
}



  def sin(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.sin(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }

  def cos(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.cos(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



def tan(v: Array[Array[Double]]): Array[Array[Double]] = {
      var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.tan(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }





def asin(v: Array[Array[Double]]): Array[Array[Double]] = {
  var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.asin(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



 def acos(v: Array[Array[Double]]): Array[Array[Double]] = {
   var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.acos(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


def   atan(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.atan(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }

def sinh(v: Array[Array[Double]]): Array[Array[Double]] = {
      var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.sinh(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



def cosh(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.cosh(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


def tanh(v: Array[Array[Double]]): Array[Array[Double]] = {
       var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
             j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.tanh(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



def exp(v: Array[Array[Double]]): Array[Array[Double]] = {
      var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.exp(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


   def log(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.log(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



def log2(v: Array[Array[Double]]): Array[Array[Double]] = {
        var conv = java.lang.Math.log(2.0)
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.log(v(i)(j))/conv
              j += 1
            }
            i += 1
        }
          vd
    }

  


def log10(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.log10(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


def  abs(v: Array[Array[Double]]): Array[Array[Double]] = {
                var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.abs(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


 def floor(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.floor(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


 def ceil(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.ceil(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }

def  transpose(Md: Array[Array[Double]]): Array[Array[Double]] =  {
    return  DoubleArray.transpose( Md ) ;
	}

  def  T(Md: Array[Array[Double]]): Array[Array[Double]] =  {
    return  DoubleArray.transpose( Md )
}



def round(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.round(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



  def sqrt(v: Array[Array[Double]]): Array[Array[Double]] = {
                var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.sqrt(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }





  def toDegrees(v: Array[Array[Double]]): Array[Array[Double]] = {
                var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.toDegrees(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



 def toRadians(v: Array[Array[Double]]): Array[Array[Double]] = {
                var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.toRadians(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


  def pow(v: Array[Array[Double]], exponent: Double): Array[Array[Double]] = {
                var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N, M)
        
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.pow(v(i)(j), exponent)
              j += 1
            }
            i += 1
        }
          vd
    }


def inv(a: Array[Array[Double]]) =  {
  scalaSci.math.LinearAlgebra.LinearAlgebra.inverse(a)
}

def eig(ar:  Array[Array[Double]])= {
    var N = ar.length
     //  the arrays of the NUMAL library are one-indexed, thus copy properly!
    var arcp = new Array[Array[Double]](N+1, N+1)
      // make a copy of the matrix to pass to comeig1
    for (r<-1 to N)
      for (c<-1 to N) 
         arcp(r)(c) = ar(r-1)(c-1)

    var  em =  new Array[Double](10)
    em(0) = 5.0e-6;   // the machine precision
    em(2)=1.0e-5;   // the relative tolerance used for the QR iteration (em(2) > em(0))
    em(4)=10.0;   // the maximum allowed number of iterations
    em(6)=1.0e-5; // the tolerance used for the eigenvectors
    em(8)=5; // the maximum allowed number of inverse iterations for the calculation of each calculation
      // arrays for the real and imaginary parts of the calculated eigenvalues of the given matrix
    var re = new Array[Double](N+1); 
    var im = new Array[Double](N+1);
      // the calculated eigenvectors are delivered in the columns of vec; an eigenvector corresponding to a real eigenvalue given in array re
      // is delivered in the corresponding column of array vec;
      // the real and imaginary part of an eigenvector corresponding to the first member of a nonreal complex conjugate pair
     // of eigenvalues given in the arrays re, im are delivered in the two consecutive columns of array vec corresponding to this pair
     // (the eigenvectors corresponding to the second members of nonreal complex complex conjugate pairs are not delivered,
     // since they are simply the complex conjugate of those corresponding to the first member of such pairs)
    var vec = new Array[Array[Double]](N+1,N+1); 
    numal.Linear_algebra.comeig1(arcp, N, em, re, im, vec )
    var vecEig= new Array[Array[Double]](N, N)  // the eigenvectors
        // the calculated eigenvectors are delivered in the columns of vec
        // one-indexed is also the [][] vec result !!
       var r = 0; var c  = 0
     while (r< N) {
       c = 0
       while  (c< N) {
          vecEig(r)(c) = vec(r+1)(c+1)  // copy eigenvectors
          c += 1
       }
       r += 1
     }

    var vecEiv= new Array[Array[Double]](N, 2)  // the eigenvalues
    r = 0
    while  (r<N)  {
        vecEiv(r)(0) = re(r+1)  // real eigenvalues
        vecEiv(r)(1) = im(r+1)  // imaginary eigenvalues
        r += 1
      }
 (vecEig, vecEiv)
}


def CholeskyL(Md: Array[Array[Double]]): Array[Array[Double]] = {
 var choleskDec = LinearAlgebra.cholesky(Md)
 var choleskMatDoubles = choleskDec.getL().getArray()
 choleskMatDoubles
}

def  Cholesky_SPD(Md: Array[Array[Double]]): Array[Array[Double]] = {
      var CholeskyDecomposition  = LinearAlgebra.cholesky(Md)
      var  choleskMatDoubles = CholeskyDecomposition.getL().getArray()
      choleskMatDoubles
}


def   det(Ma: Array[Array[Double]]) = {
     LinearAlgebra.det(Ma)
 }

/*
var tstArr:Array[Array[Double]] = Array(Array(4.0, 2.0, -1), Array(-2, 3, 7), Array(2, 3, 1))
var ev = eig(tstArr)
var mtstArr = new Matrix(tstArr)
var mev = eig(mtstArr)
*/
/*
 tstArr = [4, 2, -1;-2, 3, 7;2, 3, 1]
 [ v d ] = eig(tstArr)
*/



// computes the SVD of a given 1-indexed array  A
def svdRank(A: Array[Array[Double]]) = {
  var M = A.length
  var N = A(0).length
  var singValues = new Array[Double](N)    // the singular values
  var U = new Array[Array[Double]](M, N)
  var V = new Array[Array[Double]](N, N)

// copy data matrix to U
for (r<-1 until M)
  for (c<-1 until N)
    U(r)(c) = A(r)(c)

var em = new Array[Double](8)
em(0) = 1.0e-6  // machine's precision
em(2) = 1.0e-5  // the relative precision in the singular values
em(4) = 25.0  // the maximal number of iterations to be performed
em(5) = 1.0e-5   // the minimal nonneglectable singular value

// call the NUMAL routine to perform the SVD
 numal.Linear_algebra.qrisngvaldec(U, M-1, N-1, singValues, V, em)

// Return
// em(1): the infinity norm of the matrix
// em(3): the maximal neglected superdiagonal element
// em(5): the number of iterations performed
// em(7): the numerical rank of the matrix; i.e. the number of singular values greater than or equal to em(6)

(U,  singValues, V, em(7), em(1), em(3), em(5))
}

// computes the SVD of a given matrix A
def svd(A: Array[Array[Double]]) = {
  
var S  = svdRank(A)
  (S._1, S._2, S._3)
}



def norm2(M: Array[Array[Double]]): Double = {
  LinearAlgebra.norm2( M)
}

def  norm1(M: Array[Array[Double]]): Double = {
  LinearAlgebra.norm1( M)
}

def  normF(M: Array[Array[Double]]): Double = {
  LinearAlgebra.normF( M)
}

def  normInf(M: Array[Array[Double]]): Double = {
   LinearAlgebra.normInf( M)
}


def diag(a: Array[Double]): Array[Array[Double]] = {
         var n = a.length
         var om = new Array[Array[Double]](n, n)
         var i=1;
         while  (i<= n) {
               om(i)(i) = a(i)
               i += 1
              }
          om
    }

def  cov(v1: Array[Double], v2: Array[Double]): Double =   {
     var r = scalaSci.math.array.StatisticSample.covariance(v1, v2);
     r
 }


def covariance(v1: Array[Double], v2: Array[Double]): Double = {
  var r = cov(v1, v2)
  r
}

def  cov(v1: Array[Array[Double]], v2: Array[Array[Double]]): Array[Array[Double]] =   {
     var r = scalaSci.math.array.StatisticSample.covariance(v1, v2);
     r
 }

def  covariance(v1: Array[Array[Double]], v2: Array[Array[Double]]): Array[Array[Double]] =   {
     var r = cov(v1, v2)
     r
 }

def  cov(v: Array[Array[Double]]): Array[Array[Double]] =   {
     var r = scalaSci.math.array.StatisticSample.covariance(v, v);
     r
 }

 def  covariance(v: Array[Array[Double]]): Array[Array[Double]] =   {
     var r = scalaSci.math.array.StatisticSample.covariance(v, v);
     r
 }

def  corr(v: Array[Array[Double]]): Array[Array[Double]] =   {
     var r = scalaSci.math.array.StatisticSample.correlation(v, v);
     r
 }


def  correlation(v: Array[Array[Double]]): Array[Array[Double]] =   {
     var r =  corr(v,v);
     r
 }


def  corr(v1: Array[Double], v2: Array[Double]): Double =   {
     var r = scalaSci.math.array.StatisticSample.correlation(v1, v2);
     r
 }


def correlation(v1: Array[Double], v2: Array[Double]): Double = {
  var r = corr(v1, v2)
  r
}

def correlation(v1: Array[Array[Double]], v2: Array[Array[Double]]): Array[Array[Double]] = {
  var r = corr(v1, v2)
  r
}



def  corr(v1: Array[Array[Double]], v2: Array[Array[Double]]): Array[Array[Double]] =   {
     var r = scalaSci.math.array.StatisticSample.correlation(v1, v2);
     r
 }


def  correlation2(v1: Array[Array[Double]], v2: Array[Array[Double]]): Array[Array[Double]] = {
    var Varv1 = scalaSci.math.array.StatisticSample.variance(v1)
    var Varv2 = scalaSci.math.array.StatisticSample.variance(v2)
    var cov = scalaSci.math.array.StatisticSample.covariance(v1, v2)
    cov
}



}

object DD  {

/*
    // fast multiplication for 2-dimensional arrays
  def *(A: Array[Array[Double]], B: Array[Array[Double]]): Array[Array[Double]] = {
    var rN= A.length;  var rM = A(0).length
    var  sN = B.length;  var sM = B(0).length

   var  v1Colj = new Array[Double](rM)
   var result = new Array[Array[Double]](rN, sM)
   var j=0; var k=0;
   while (j < sM)  {
       k=0
      while  (k < rM) {
        v1Colj(k) =B(k)(j)
        k += 1
      }

      var i=0;
      while (i<rN) {
        var   Arowi = A(i)
        var   s = 0.0;
        k=0
        while (k< rM) {
          s += Arowi(k)*v1Colj(k)
          k += 1
        }
      result(i)(j) = s;
      i += 1
      }
 j += 1
   }
  result
}
*/

 def *(A: Array[Array[Double]], B: Array[Array[Double]]): Array[Array[Double]] = {
    var rN= A.length;  var rM = A(0).length
    var  sN = B.length;  var sM = B(0).length

   var  v1Colj = new Array[Double](rM)
   var result = new Array[Array[Double]](rN, sM)
   for (j<-0 to sM-1)  {
      for (k<-0 to rM-1) 
        v1Colj(k) =B(k)(j)

      for ( i<-0 to rN-1)  {
        var   Arowi = A(i)
        var   s = 0.0;
        for (k<-0 to rM-1) 
          s += Arowi(k)*v1Colj(k)
        
      result(i)(j) = s;
      }
   }

  result
}



// fast addition for 2-dimensional arrays
  def +(A: Array[Array[Double]], B: Array[Array[Double]]): Array[Array[Double]] = {
    var  rN= A.length;  var rM = A(0).length
    var  sN = B.length;  var sM = B(0).length

    if (rN != sN || rM != sM) {
      println("dimensions of matrices are not equal in matrix addition")
      return null
    }
    var result = new Array[Array[Double]](rN, rM)
    for  (r<-0 to rN-1 )
      for  (c<-0 to rM-1) 
        result(r)(c) = A(r)(c)+B(r)(c)
       
    return result
  }


// fast subtraction for 2-dimensional arrays
  def -(A: Array[Array[Double]], B: Array[Array[Double]]): Array[Array[Double]] = {
    var  rN= A.length;  var rM = A(0).length
    var  sN = B.length;  var sM = B(0).length

    if (rN != sN || rM != sM) {
      println("dimensions of matrices are not equal in matrix subtraction")
      return null
    }
    var result = new Array[Array[Double]](rN, rM)
    for (r<-0 to rN-1) 
       for (c<-0 to rM-1) 
        result(r)(c) = A(r)(c)-B(r)(c)
    return result
  }

def *(x: Double, A: Array[Array[Double]]): Array[Array[Double]] = {
    var  N= A.length;  var M = A(0).length

    var result = new Array[Array[Double]](N, M)
    for  (r<-0 to N-1 ) 
       for (c<-0 to  M-1 ) 
        result(r)(c) = A(r)(c)*x
    return result
  }


def +(x: Double, A: Array[Array[Double]]): Array[Array[Double]] = {
    var  N= A.length;  var M = A(0).length

    var result = new Array[Array[Double]](N, M)
    for (r<-0 to N-1) 
      for (c<-0 to M-1)  
        result(r)(c) = A(r)(c)+x
        
    return result
  }


def -(x: Double, A: Array[Array[Double]]): Array[Array[Double]] = {
    var  N= A.length;  var M = A(0).length

    var result = new Array[Array[Double]](N, M)
    for (r<-0 to N-1 ) 
      for (c<-0 to M-1 ) 
        result(r)(c) = A(r)(c)-x
    return result
  }

def *(A: Array[Array[Double]],  x: Double):Array[Array[Double]] =  { 
    var  N= A.length;  var M = A(0).length

    var result = new Array[Array[Double]](N, M)
    for  (r<-0 to N-1 ) 
       for (c<-0 to  M-1 ) 
        result(r)(c) = A(r)(c)*x
    return result
}

def +(A: Array[Array[Double]],  x: Double):Array[Array[Double]] =  { 
  var  N= A.length;  var M = A(0).length

    var result = new Array[Array[Double]](N, M)
    for (r<-0 to N-1) 
      for (c<-0 to M-1)  
        result(r)(c) = A(r)(c)+x
        
    return result
  }


def -(A: Array[Array[Double]],  x: Double):Array[Array[Double]] =  { 
    var  N= A.length;  var M = A(0).length

    var result = new Array[Array[Double]](N, M)
    for (r<-0 to N-1 ) 
      for (c<-0 to M-1 ) 
        result(r)(c) = A(r)(c)-x
    return result
  }




  def randf(n: Int, m:Int): Array[Array[Double]] = {
      var  v =new Array[Array[Double]](n,m)
      for (i<-0 to  n-1) 
           for (j <-0 to  m-1) 
              v(i)(j) = java.lang.Math.random
        v
    }




}
